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Long ranged electrostatic interactions are time consuming to calculate in molecular dynamics and Monte- 
Carlo simulations. We introduce an algorithmic framework for simulating charged particles which modifies 
the dynamics so as to allow equilibration using a JocaJ Hamiltonian. The method introduces an auxiliary field 
with constrained dynamics so that the equilibrium distribution is determined by the Coulomb interaction. We 
demonstrate the efficiency of the method by simulating a simple, charged lattice gas. 
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The electrostatic interaction between two point charges 
in a medium with uniform dielectric constant eq varies as 
eie2/4:TTeQr. The large numerical value of this energy together 
with its long range are such that it is very often the most costly 
component in the simulation of charged condensed matter sys- 
tems. Naive evaluation of the electrostatic energies in molecu- 
lar dynamics and Monte-Carlo algorithms leads to inner loops 
where the summation over all pairs takes a time which scales 
as 0{N'^) for a single step in which all N particles are up- 
dated. 

Many methods are used to improve this poor scaling: The 
optimized Ewald algorithm splits the summation between real 
and Fourier space and has a complexity of 0{N'^/'^) \ l,|2j]. By 
interpolation of the charge distribution onto a grid, fast Fourier 
transform methods allow a scaling in 0(7Vlog(A^)) 0]. Fi- 
nally, a popular method in very large simulations is an expan- 
sion of the charge distribution using hierarchical multipoles 
|4, 5]. The asymptotic improvement in efficiency comes, how- 
ever, with great increases in the complexity of the coding, es- 
pecially when distributed on multiprocessor computers. The 
numerical prefactors in these scaling laws are uncomfortably 
high: Despite the great effort put into optimizing the elec- 
trostatic loop, it is found that in the simulation of a large 
biomolecule (with N ^ 10^) the great majority of the CPU 
time is still used in the Coulomb loop |6] in even the most so- 
phisticated numerical codes. Most of these "fast" methods can 
only be used efficiently in molecular dynamics simulations; 
there are many occasions where one would like to perform ef- 
ficient Monte-Carlo simulations due to the stability and sim- 
plicity of the method. 

The classical methods for treating charged systems have an- 
other disadvantage, their inability to treat systems with inho- 
mogeneous dielectric constants. Dielectric inhomogeneities 
have drastic effects on material properties. For instance the di- 
electric contrast between water and the core of proteins leads 
to expulsion of counter-ions from a thick hydration layer 
0]. To treat charging effects in proteins quite arbitrary, uncon- 
trolled approximations are made |8] on effective electrostatic 
interactions in the vicinity of a protein in order to reduce in- 
teractions to effective pair-wise additive potentials. Similarly 
much work |9| has been performed on the phase structure of 
charged synthetic polymers while neglecting the large dielec- 
tric contrasts between water based solvents and oily backbone 



structures which are surely important in the discussion of the 
stability of the necklace structures predicted in these systems. 
At present the most promising algorithms are based on the 
non-local Marcus energy functional L10i..l 1.1 . 

This letter introduces a local algorithm with a propagating 
field E with purely local dynamics on an interpolating grid; 
it has a complexity in 0{N) and is elementary to implement. 
In contrast to conventional grid methods we do not solve for 
all the field variables at each integration step; we let the field 
evolve with its own intrinsic dynamics. We were motivated 
by the observation that Maxwell's equations, which are lo- 
cal, produce Coulomb interactions due to the propagation of a 
vectorial field. These dynamic equations are, as we shall see, 
not the only dynamic way of generating the Coulomb inter- 
action. Our method allows a direct, local implementation of 
dielectric inhomogeneities. As a demonstration of the method 
we present an explicit implementation of a local Monte-Carlo 
algorithm for a charged lattice gas. We note that techniques 
which interpolate charge degrees of freedom onto a lattice are 
already very well understood; they form part of standard pack- 
ages such as Amber |3]. 

We proceed by showing that the Coulomb interaction can 
be derived from a constrained variational problem. We then 
show that the constraint equations are solved locally if we al- 
low electric fields which have both gradient and rotational de- 
grees of freedom. This freedom can be used to produced a 
local Monte-Carlo algorithm. Finally we present a numerical 
verification of the method. 

The energy of a system of charged particles in a uniform 
dielectric background is expressed as a function of the electric 
field E 

U^eoJ^d^r (1) 

where the electric field is constrained by Gauss's law 

div E - p/ea = (2) 

It is known from classical electrostatics that one solution of 
equation is given by E = — grad so that = —p/eo- 
The general solution to the constraint eq. (|2} is thus 

E = -grad </« + curl Q (3) 
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where is unique to within an additive constant and Q arbi- 
trary. In Fourier space the electric field can be expressed as 
E(k) = — ik0 + zk A Q. The second term of this expression 
is perpendicular to k so that there are two physical degrees 
of freedom in the Q field, corresponding to two independent 
polarization states. We can consider that the field is due to a 
static potential plus transverse photons. 

Let us study the stationary states to the variational problem 
posed by eq. Q and eq. Q by using a Lagrange multiplier, 
with the functional 



E 

- A(r)(eodiv E - p) 



d^r 



(4) 



implying that E + grad A = 0. The Lagrange multiplier is 
identical to the static electrostatic potential, </> and the min- 
imum energy is Ucouiomb = ^ / (grad 0)^ d'^r. Consider 
the energy eq. Q for an arbitrary E satisfying the constraints 
then 

^ = f / [(grad cj^f + (curl Q)2] d^r (5) 

Cross terms vanish, as is shown by integrating by parts. 

We now turn to the statistical mechanics of a field with the 
energy of eq. Q constrained by Gauss's law. We do not im- 
pose that the electric field is calculated from a potential. The 
partition function of a fixed set of charges in presence of a 
fluctuating field E is given by 



Z{{v})= / PEe 



n,5(divE-p(r)/eo) (6) 



The argument {r} denotes the fact that the integral over the 
particle positions has not yet been performed. There are two 
ways of treating this equation. Either one introduces an inte- 
gral representation of the delta function or one notes that the 
integral over the field E decomposes into a (unique) gradient 
term and a (non unique) rotation so that 



Z{{v}) 



I (grad ( 



d-^r 



/ d^r 



(7) 



where PEf = Hr '^('ii^ E)2?E performs the summation over 
all the rotational degrees of freedom of the field described by 
the potential Q. All the dependence on the particle positions 
is in the prefactor characterized by the electrostatic potential 
(f) found by solving Poisson's equation. This prefactor gives 
the Coulomb interaction between the particles. The remaining 
integral is independent of the positions of the charges; rather 
remarkably integration over the full set of fields allowed by 
the constraint multiplies the standard partition function by a 
simple constant. This extra factor in the partition function can 
be ignored. 

In the presence of non-uniform dielectric media one pro- 
ceeds in a similar manner with the energy 



26(r) 



d^r 



(8) 



and the constraint div D p = 0. We deduce that the 
displacement is given by D = — egrad^ + curl Q with 
div (egrad (/>) = — p and 



so that 



-f / e(r)(grad 4>f d'r J ^-^^ ^- J d^ 



(9) 



Z({r}) = Zcoulornbi{r})Zfiuct{{r}) 



(10) 



This time the normalization is a function of the distribution of 
dielectric inhomogeneities. When implemented in inhomoge- 
neous media our treatment leads to potentials which are the 
sum of the Coulomb and a fluctuation potential |12] which 
varies as l/r^ for two widely separated particles. This term 
comes from thermally driven dipole-dipole interactions: Fluc- 
tuations in the field produce an inhomogeneous polarization, 
P, of the dielectric background. This produces an equivalent 
charge density of —div P which interacts via Coulomb's law. 
Such fluctuation potentials are to be expected from the Lif- 
shitz theory of dielectrics. 

We now propose a lattice version of the above equa- 
tions suitable for numerically studying the thermodynamics 
of charged systems. The trick is to use the arbitrary vector po- 
tential Q to simplify the calculation of the updated fields after 
the motion of a charged particle. We need, also, to sum over 
all rotational degrees of freedom of the field in order to calcu- 
late statistical weights from the partition function eq. (|6j. 
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FIG. 1 : Top: A charged particle is present on the left lattice point. If 
the particle is transferee! to the right lattice point the constraints are 
still satisfied if £1,2 is modified to £1,2 — e/eo on the connecting 
link where e is the charge of the particle. Bottom: The four fields 
associated with a single plaquette, C, are modified by a rotational 
motion. 

The observation that we use in order to implement a local 
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algorithm is that the constraint equation, eq. (|2j can be up- 
dated locally in a system in which charge is conserved. We 
firstly reinterpret the constraint in terms of Faraday's concept 
of conserved electric flux: Consider, fig. H (top), a network 
where the charges are confined to the vertices {i} and the field 
is associated with links between two sites, Around each 

lattice point we imagine a cube and write the constraint in in- 
tegral form j E.dS = Ci/eo. The integral is over the surface 
of the cube, is the enclosed charged at the site. We use 
the notation i?i 2 to denote the total flux leaving 1 towards 2; 
clearly £'12 = — i?2,i. 

The discretized version of the integral constraint is 
J2j = Ci/eQ The discretized energy is given by 

links 

Here and in what follows we assume that the lattice spacing is 
unity. 

Start with a system where the constraint is satisfied, fig.H 
(top), and displace a charge, e, situated on the leftmost lattice 
site, 1, to the rightmost site, 2. The constraint is again satisfied 
at both sites if the field associated with the connecting link is 
updated according to the rule i?i 2 i?i,2 ^e/^o- This is our 
Monte-Carlo move for the particles, involving a correlated up- 
date of a single charge and the field on the link connecting two 
sites. To update the field configurations fig.[n(bottom) we up- 
date all the field values of a single plaquette while conserving 
the constraint at each vertex. In fig. [2 (bottom) £'1.2 and £'4.1 
increase by an increment A whereas £4 3 and £3 2 decrease 
by A so that at each vertex the sum of the entering and leav- 
ing fields is constant. It is this last update that performs the 
integration over all the rotational degrees of freedom in the E 
field. 

The two moves are not quite sufficient to equilibrate a sys- 
tem with periodic boundary conditions in all situations. This 
problem is linked with the solution (f> = — E.r or E = E of 
the Laplace equation on a torus where E is an arbitrary con- 
stant vector Motion of the charges generates fluctuations in 
E, while updates such as those in fig.n(bottom) preserve E; 
similar phenomena also occur with Maxwell's equations. In 
order to be absolutely sure that the algorithm is ergodic we 
introduced a third possible Monte-Carlo step which consists 
of a shift in E. By keeping track of the evolution of E as the 
particles move this last update can be efficiently implemented 
without destroying the 0{N) scaling of the algorithm. In the 
largest systems fluctuations in this single mode should give a 
small contribution to the thermodynamics if the initial condi- 
tion is typical; in such cases this update can be eliminated. 

We have performed two initial verifications of the algo- 
rithm. Firstly we randomly placed four positive and four neg- 
ative charged particles on a 4 x 4 x 4 lattice with periodic 
boundary conditions. We performed field updates using the 
Metropolis algorithm at zero temperature in order to anneal 
the field E. We then solved for the electric fields using a 
standard linear algebra package. The results were identical to 



within numerical errors. Annealing E was crucial in order to 
get agreement between the two methods with frozen charges. 
A second check was then performed with two different im- 
plementations of the Metropolis algorithm now run at several 
finite temperatures: 36 mutually avoiding, charged particles 
were distributed in a 6 x 6 x 6 cube with periodic boundary 
conditions. The dielectric constant was uniform. In the first 
simulation the linear solver was used to calculate the exact in- 
teraction energy of the charged particles at each Monte-Carlo 
step. In the second simulation we implemented the above lo- 
cal algorithm. For the two methods we then compared the 
static structure factor for charge-charge correlations, finding 
good agreement. 

We verified the efficiency of the algorithm by determin- 
ing the autocorrelation time of the slowest density and charge 
modes in a cube of dimension L using the method described 
in fig. (EJ. We plot the relaxation time in Monte-Carlo 
sweeps; during each sweep an average of one Monte-Carlo 
trial is performed for each degree of freedom in the simula- 
tion. As expected in a Monte-Carlo algorithm dominated by 
diffusive motion the slowest mode relaxes in a time which 
varies as L^. There is no anomalous or critical slowing down 
due to the coupling between the particles and the electric field. 
The saturation of the relaxation time for charge fluctuations at 
large £ is a consequence of screening; charge fluctuations re- 
lax by diffusion over the Debye length. 
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FIG. 2: Relaxation times in Monte-Caiio sweeps for the density- 
density (top) correlation function and the charge-charge (bottom) 
correlation function measured for the mode ^(1,0, 0) plotted as a 
function of L^. The top curve shows scaling compatible with simple 
diffusion. The bottom curve saturates for large systems sizes. The 
multiple points for each L correspond to three independent determi- 
nations of the relaxation time. System sizes between L = 6 and 
L — 36 with charged particles. 6 x 10^ sweeps per simulation 
for a total simulation time of three days on a AMD Athlon computer. 

We have shown that the thermodynamics of charged sys- 
tems can be simulated locally by introducing a propagating 
vector field so that particles interact via retarded, diffusing 
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fields. The dynamic properties of the system are strongly 
modified but by construction thermodynamics is an invariant 
of the propagation dynamics. Our treatment of the dynamics 
of the field E is similar to the Coulomb or radiation gauge in 
classical electrodynamics: Normally one writes that 



dA 

E = -grad0- — 



(12) 



where in the Coulomb gauge div A 0. There is a rather 
close analogy between our dynamic scheme of solving for the 
constraint equations with certain methods of quantization in 
the Coulomb gauge: In Dirac's quantization of the electro- 
dynamic field Gauss's law is a weak identity fl4ll dependent 
on the choice of the initial wave function. In our simulation 
Gauss's law is the result of a restricted choice of possible 
moves in the Monte-Carlo algorithm together with a special 
initial condition. 

What are the advantages of the present method over direct 
integration of Maxwell's equations which are also an example 
of an 0{N) algorithm [18]? Monte-Carlo algorithms are par- 
ticularly easy to implement and have good stability with large 
step sizes. In addition, we have checked that the fast equilibra- 
tion of the electric degrees of freedom in fig. (2) allows one 
to perform field updates far more rarely than particle updates 
leading to additional acceleration of the algorithm. Such mul- 
tiple time step ideas have been applied to conventional electro- 
static solvers but in molecular dynamics are sometimes prone 
to numerical instabilities jQ]. 

In the implementation of the algorithm we were inspired by 
recent work on hydrodynamic interactions (varying as l/r) 
via a Lattice-Boltzmann algorithm in simulations of polymer 
solutions 1 15]. Another analogous problem to ours has been 
treated by Car and Parrinello who have shown 1 16] that intro- 
duction of a fictitious dynamics leads to much improved effi- 
ciencies in solving for constraints; our work differs in that the 
constraint of Gauss's law is solved exactly at each simulation 
step whereas the Car-Parrinello algorithm leads to a simulated 
annealing solution of the constraint equations Q. 

Finally let us note that other functionals do exist for the 
electric potential in the presence of sources iITtIi . In particular 
the functional 



seems, at first sight particularly simple and attractive. Unfor- 
tunately, the minimum of this functional is minus the correct 
electrostatic energy. It can not be used as a functional for both 
the field evolution and the particle motion. Application of the 
algorithm to large atomistic systems remain to be tested, but 
our method provides an alternative to existing treatments of 
Coulomb interactions. 
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